TFG def

Librerias y paquetes

library(zoo)

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(xts)
library(CASdatasets)

## Loading required package: survival

library(dplyr)

## 
## ######################### Warning from 'xts' package ##########################
## #                                                                             #
## # The dplyr lag() function breaks how base R's lag() function is supposed to  #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
## # source() into this session won't work correctly.                            #
## #                                                                             #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
## # dplyr from breaking base R's lag() function.                                #
## #                                                                             #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## ###############################################################################

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:xts':
## 
##     first, last

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(tidyr)
library(ggplot2)
library(scales)
library(forcats)
library(kableExtra)

## 
## Attaching package: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

library(survival)
library(effects)

## Loading required package: carData

## lattice theme set by effectsTheme()
## See ?effectsTheme for details.

library(MASS)
library(broom)

## Warning: package 'broom' was built under R version 4.5.2

library(ggplot2)
library(mgcv)

## Loading required package: nlme

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:dplyr':
## 
##     collapse

## This is mgcv 1.9-4. For overview type '?mgcv'.

library(stringr)
library(statmod) 
library(tweedie)
library(openxlsx)
library(glmmTMB)

## Warning in check_dep_version(dep_pkg = "TMB"): package version mismatch: 
## glmmTMB was built with TMB package version 1.9.17
## Current TMB package version is 1.9.19
## Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)

## 
## Attaching package: 'glmmTMB'

## The following object is masked from 'package:statmod':
## 
##     tweedie

library(cplm)

## Loading required package: coda

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

## Loading required package: splines

## 
## Attaching package: 'cplm'

## The following object is masked from 'package:glmmTMB':
## 
##     VarCorr

## The following object is masked from 'package:nlme':
## 
##     VarCorr

Instalación y carga de datos

#install.packages('xts')
install.packages("CASdatasets", repos = "http://dutangc.free.fr/pub/RRepos/pub/", type="source")

Carga de los 2 datasets

# Tabla frecuencia
data(freMTPL2freq)
freq <- freMTPL2freq

# Tabla severidad
data(freMTPL2sev)
sev <- freMTPL2sev

Exploración inicial de la estructura de la base de datos

Tabla de frecuencia

Exploración inicial –>

summary(freq)

##      IDpol           ClaimNb          Exposure           VehPower     
##  1      :     1   Min.   : 0.000   Min.   :0.002732   Min.   : 4.000  
##  3      :     1   1st Qu.: 0.000   1st Qu.:0.180000   1st Qu.: 5.000  
##  5      :     1   Median : 0.000   Median :0.490000   Median : 6.000  
##  10     :     1   Mean   : 0.039   Mean   :0.528743   Mean   : 6.455  
##  11     :     1   3rd Qu.: 0.000   3rd Qu.:0.990000   3rd Qu.: 7.000  
##  13     :     1   Max.   :16.000   Max.   :2.010000   Max.   :15.000  
##  (Other):677985                                                       
##      VehAge           DrivAge        BonusMalus        VehBrand     
##  Min.   :  0.000   Min.   : 18.0   Min.   : 50.00   B12    :166022  
##  1st Qu.:  2.000   1st Qu.: 34.0   1st Qu.: 50.00   B1     :162730  
##  Median :  6.000   Median : 44.0   Median : 50.00   B2     :159854  
##  Mean   :  7.044   Mean   : 45.5   Mean   : 59.76   B3     : 53393  
##  3rd Qu.: 11.000   3rd Qu.: 55.0   3rd Qu.: 64.00   B5     : 34752  
##  Max.   :100.000   Max.   :100.0   Max.   :230.00   B6     : 28545  
##                                                     (Other): 72695  
##      VehGas       Area          Density     
##  Diesel :332128   A:103952   Min.   :    1  
##  Regular:345863   B: 75457   1st Qu.:   92  
##                   C:191874   Median :  393  
##                   D:151592   Mean   : 1792  
##                   E:137163   3rd Qu.: 1658  
##                   F: 17953   Max.   :27000  
##                                             
##                          Region      
##  Centre                     :160595  
##  Rhone-Alpes                : 84748  
##  Provence-Alpes-Cotes-D'Azur: 79315  
##  Ile-de-France              : 69789  
##  Bretagne                   : 42120  
##  Pays-de-la-Loire           : 38747  
##  (Other)                    :202677

names(freq)

##  [1] "IDpol"      "ClaimNb"    "Exposure"   "VehPower"   "VehAge"    
##  [6] "DrivAge"    "BonusMalus" "VehBrand"   "VehGas"     "Area"      
## [11] "Density"    "Region"

dim(freq)

## [1] 677991     12

12 variaable y 677991 registros.

Renombramiento variables –>

freq <- freq %>%
  rename(
    id_poliza = IDpol,
    exposicion = Exposure,
    numero_siniestros = ClaimNb,
    potencia = VehPower,
    edad_conductor = DrivAge,
    edad_vehiculo = VehAge,
    bonus_malus = BonusMalus,
    marca = VehBrand,
    combustible = VehGas,
    area = Area,
    region = Region
  )

Tabla de severidad

Exploración inicial –>

summary(sev)

##      IDpol        ClaimAmount     
##  2241683:   16   Min.   :      1  
##  3253234:   11   1st Qu.:    686  
##  3254353:   11   Median :   1172  
##  2248174:    9   Mean   :   2266  
##  2239279:    8   3rd Qu.:   1212  
##  2216294:    6   Max.   :4075401  
##  (Other):26383

names(sev)

## [1] "IDpol"       "ClaimAmount"

dim(sev)

## [1] 26444     2

2 variables y 26444 registros

Renombramiento variables –>

sev <- sev %>%
  rename(
    id_poliza = IDpol,
    coste_siniestro = ClaimAmount
  )

Verificación básica de la calidad de los datos

colSums(is.na(freq))

##         id_poliza numero_siniestros        exposicion          potencia 
##                 0                 0                 0                 0 
##     edad_vehiculo    edad_conductor       bonus_malus             marca 
##                 0                 0                 0                 0 
##       combustible              area           Density            region 
##                 0                 0                 0                 0

colSums(is.na(sev))

##       id_poliza coste_siniestro 
##               0               0

Analisis valores duplicados –>

any(duplicated(freq$IDpol))   

## [1] FALSE

any(duplicated(sev$IDpol))   

## [1] FALSE

Preparación y unificación del conjunto de datos

sev <- sev %>%
  group_by(id_poliza) %>%
  summarise(
    total_importe_siniestros = sum(coste_siniestro, na.rm = TRUE)
  ) %>%
  ungroup()   

Unificamos las dos tablas

polizas <- freq %>%
  left_join(sev, by = "id_poliza")

Remplazamos los NA en el importe total de siniestros y en el número de siniestros por cero

polizas <- polizas %>%
  mutate(
    total_importe_siniestros = ifelse(is.na(total_importe_siniestros), 0, total_importe_siniestros),
  )

Variables suplementarias

Tasa de siniestralidad

Cuántos siniestros ocurren en promedio por póliza, ajustado por el tiempo que la póliza ha estado en vigor. Permite comprar pólizas que no han estado aseguradas el mismo tiempo.

\[ \text{tasa siniestros} = \frac{\text{numero siniestros}}{\text{exposición}} \]

polizas <- polizas %>%
  mutate(tasa_siniestros = numero_siniestros / exposicion)

Prima pura estimada

Cantidad de dinero que, en promedio, se espera pagar en siniestros por póliza y unidad de exposición.

\[ \text{importe medio siniestro} = \begin{cases} \dfrac{\text{total importe siniestros}}{\text{numero siniestros}}, & \text{si } \text{numero siniestros} > 0 \\[6pt] 0, & \text{si } \text{numero siniestros} = 0 \end{cases} \]

\[ \text{prima pura estimada} = \text{tasa siniestros} \times \text{importe medio siniestro} \]

polizas <- polizas %>%
  mutate(
    importe_medio_siniestro = ifelse(numero_siniestros > 0, total_importe_siniestros / numero_siniestros, 0),
    prima_pura_estimada = tasa_siniestros * importe_medio_siniestro
  )

Analisis descriptivo

Exposición

summary(polizas$exposicion)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.002732 0.180000 0.490000 0.528743 0.990000 2.010000

Distribución –>

ggplot(polizas, aes(x = exposicion)) +
  geom_histogram(bins = 10, fill = "steelblue", color = "white") +
  labs(title = "Distribución de la exposición", x = "Años asegurados", y = "Número de pólizas")+
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 12),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8)
  )

ggplot(polizas, aes(x = "", y = exposicion)) +
  geom_boxplot(fill = "steelblue", color = "grey2") +
  labs(
    title = "Distribución de la exposición",
    y = "Exposición (años asegurados)",
    x = ""
  ) +
  theme_minimal()+
  theme(
    plot.title = element_text(hjust = 0.5, size = 12),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8)
  )

Depuramos los valores superiores a un año –>

# Modificar las observaciones con exposición > 1, igualándolas a 1
polizas <- polizas %>%
  mutate(exposicion = ifelse(exposicion > 1, 1, exposicion))

# Verificar el cambio
summary(polizas$exposicion)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.002732 0.180000 0.490000 0.528537 0.990000 1.000000

Distribución después de depurar la variable

ggplot(polizas, aes(x = exposicion)) +
  geom_histogram(bins = 10, fill = "steelblue", color = "white") +
  labs(title = "Distribución de la exposición", x = "Años asegurados", y = "Número de pólizas")+
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 12),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8)
  )

Número de siniestros

Tabla descriptiva antes de la depuración

intervalos <- c(-Inf, 0, 1, 2, 3, 5, 10, 20, Inf)
etiquetas  <- c("0", "1", "2", "3", "4-5", "6-10", "11-20", "Más de 20")

tabla_siniestros <- table(cut(polizas$numero_siniestros, breaks = intervalos, labels = etiquetas, right = TRUE)) %>%
  as.data.frame()

colnames(tabla_siniestros) <- c("Intervalo", "Recuento")
tabla_siniestros$Porcentaje <- round(tabla_siniestros$Recuento / sum(tabla_siniestros$Recuento) * 100, 2)

tabla_siniestros

##   Intervalo Recuento Porcentaje
## 1         0   653047      96.32
## 2         1    23571       3.48
## 3         2     1298       0.19
## 4         3       62       0.01
## 5       4-5        7       0.00
## 6      6-10        3       0.00
## 7     11-20        3       0.00
## 8 Más de 20        0       0.00

Eliminamos los valores superiores a 4 siniestros para evitar una sobre dispersión.

polizas <- polizas %>%
  mutate(numero_siniestros = ifelse(numero_siniestros > 4, 4, numero_siniestros))
summary(polizas$numero_siniestros)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.03895 0.00000 4.00000

ggplot(polizas, aes(x = factor(numero_siniestros))) +
  geom_bar(fill = "steelblue", color = "white") +
  labs(
    title = "Número de siniestros por póliza",
    x = "Número de siniestros",
    y = "Frecuencia"
  )  +
  theme_minimal(base_size = 10) +
  theme(
    plot.title    = element_text(hjust = 0.5, face = "bold", size = 12),
    plot.subtitle = element_text(hjust = 0.5, size = 9)
  )

tabla_sin <- polizas %>% 
  count(numero_siniestros) %>% 
  mutate(
    prop = n / sum(n),
    etiqueta = percent(prop, accuracy = 0.1, decimal.mark = ","),
    numero_siniestros = factor(numero_siniestros)
  )

ggplot(tabla_sin, aes(x = numero_siniestros, y = n)) +
  geom_col(fill = "steelblue", color = "white") +
  geom_text(aes(label = etiqueta), vjust = -0.3, size = 3) +
  labs(
    title    = "Número de siniestros por póliza",
    x = "Número de siniestros en la póliza",
    y = "Número de pólizas"
  ) +
  scale_y_continuous(
    labels = label_number(big.mark = ".", decimal.mark = ","),
    expand = expansion(mult = c(0, 0.05))
  ) +
  theme_minimal(base_size = 10) +
  theme(
    plot.title    = element_text(hjust = 0.5, face = "bold", size = 12),
    plot.subtitle = element_text(hjust = 0.5, size = 9)
  )

#### Tasa de sineistros

Proporción de pólizas con tasa igual a cero

tabla_cero <- polizas %>%
  mutate(tasa_pos = if_else(tasa_siniestros == 0, "tasa siniestros = 0", "tasa sineistros > 0 siniestros")) %>%
  count(tasa_pos) %>%
  mutate(porcentaje = 100 * n / sum(n))

ggplot(tabla_cero, aes(x = tasa_pos, y = n)) +
  geom_col(fill = "steelblue", color = "white") +
  geom_text(aes(label = paste0(round(porcentaje, 1), "%")),
            vjust = -0.5, size = 3) +
  labs(
    title = "Pólizas con tasa de siniestros nula y positiva",
    x = "Situación",
    y = "Número de pólizas"
  )+
  theme_minimal(base_size = 10) +
  theme(
    plot.title    = element_text(hjust = 0.5, face = "bold", size = 12),
    plot.subtitle = element_text(hjust = 0.5, size = 9)
  )

Histograma solo para tasas positivas y recortado

tasa_pos <- polizas %>%
  filter(tasa_siniestros > 0,
         tasa_siniestros <= 5)  # recorte a 5; ajusta si hace falta

ggplot(tasa_pos, aes(x = tasa_siniestros)) +
  geom_histogram(bins = 50, fill = "steelblue", color = "white") +
  labs(
    title = "Distribución de la tasa de siniestros (pólizas con tasa > 0,\nrecorte en 5)",
    x = "Tasa de siniestros",
    y = "Número de pólizas"
  ) +
  theme_minimal(base_size = 10) +
  theme(
    plot.title    = element_text(hjust = 0.5, face = "bold", size = 12),
    plot.subtitle = element_text(hjust = 0.5, size = 9)
  )

tasa_pos <- polizas %>% 
  filter(tasa_siniestros > 0,
         tasa_siniestros <= quantile(tasa_siniestros[tasa_siniestros > 0], 0.99))

ggplot(tasa_pos, aes(x = tasa_siniestros, y = exposicion)) +
  geom_point(alpha = 0.15, color = "steelblue", size = 1) +
  labs(
    title = "Dispersión de la tasa de siniestros frente a la exposición",
    x = "Tasa de siniestros",
    y = "Exposición"
  ) +
  theme_minimal(base_size = 10) +
  theme(
    plot.title    = element_text(hjust = 0.5, face = "bold", size = 12),
    plot.subtitle = element_text(hjust = 0.5, size = 9)
  )

Importe total siniestro

ggplot(polizas %>% filter(total_importe_siniestros > 0),
       aes(x = total_importe_siniestros)) +
  geom_histogram(bins = 50, fill = "steelblue", color = "white") +
  scale_x_log10(
    labels = comma_format(accuracy = 1)) +
  labs(title = "Importe total de siniestros (log)", x = "€ (escala log)", y = "Frecuencia") +
  theme_minimal(base_size = 10) +
  theme(
    plot.title    = element_text(hjust = 0.5, face = "bold", size = 12),
    plot.subtitle = element_text(hjust = 0.5, size = 9)
  )

####Importe medio de siniestros

ggplot(polizas %>% filter(importe_medio_siniestro > 0),
       aes(x = importe_medio_siniestro)) +
  geom_histogram(bins = 50, fill = "steelblue", color = "white") +
  scale_x_log10(
    labels = comma_format(accuracy = 1)) +
  labs(title = "Importe medio de siniestro (log)", x = "€ (escala log)", y = "Frecuencia") +
  theme_minimal(base_size = 10) +
  theme(
    plot.title    = element_text(hjust = 0.5, face = "bold", size = 12),
    plot.subtitle = element_text(hjust = 0.5, size = 9)
  )

####Prima pura estimada

options(scipen = 999)  # evita usar notación científica
summary(polizas$prima_pura_estimada)

##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
##        0.0        0.0        0.0      383.3        0.0 18524548.0

polizas %>% 
  filter(prima_pura_estimada > 0) %>%
  ggplot(aes(x = prima_pura_estimada)) +
  geom_histogram(bins = 60, fill = "steelblue", color = "white") +
  scale_x_log10(labels = label_number(big.mark = ".", decimal.mark = ",")) +
  labs(title = "Prima pura estimada (>0, escala log10)", x = "€ (log10)", y = "Frecuencia") +
  theme_minimal(base_size = 10) +
  theme(
    plot.title    = element_text(hjust = 0.5, face = "bold", size = 12),
    plot.subtitle = element_text(hjust = 0.5, size = 9)
  )

Edad del conductor

Distribución edad del conductor

ggplot(polizas, aes(x = edad_conductor)) +
  geom_histogram(bins = 50, fill = "steelblue", color = "white") +
  labs(title = "Distribución de la edad del conductor", x = "Edad", y = "Frecuencia")+
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 12),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8)
  )

Estudiamos cómo categorizar los datos:

exposicion_por_edad <- polizas %>%
  group_by(edad_conductor) %>%
  summarise(exposicion_total = sum(exposicion, na.rm = TRUE))

ggplot(exposicion_por_edad, aes(x = edad_conductor, y = exposicion_total)) +
  geom_col(fill = "steelblue", color = "white") +
  labs(
    title = "Distribución de la exposición según la edad del conductor",
    x = "Edad del conductor",
    y = "Exposición"
  )  +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 12, face = "bold"),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8)
  )

siniestros_por_edad <- polizas %>%
  group_by(edad_conductor) %>%
  summarise(media_siniestros = mean(numero_siniestros, na.rm = TRUE))

ggplot(siniestros_por_edad, aes(x = edad_conductor, y = media_siniestros)) +
  geom_col(fill = "steelblue", color = "white") +
  labs(
    title = "Número medio de siniestros por edad del conductor",
    x = "Edad del conductor",
    y = "Número medio de siniestros"
  )  +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 12, face = "bold"),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8)
  )

Categorizamos la variable:

polizas <- polizas %>%
  mutate(
    edad_conductor_cat = cut(
      edad_conductor,
      breaks = c(-Inf, 25, 40, 60,  Inf),
      labels = c("Jóvenes", "Adultos jóvenes", "Adultos","Mayores"),
      right = TRUE
    )
  )

Estudiamos la variable edad categorizada:

porcentajes_edad <- polizas %>%
  count(edad_conductor_cat) %>%
  mutate(porcentaje = round(100 * n / sum(n), 1))

ggplot(porcentajes_edad, aes(x = edad_conductor_cat, y = n)) +
  geom_col(fill = "steelblue", color = "white") +
  geom_text(aes(label = paste0(porcentaje, "%")), 
            vjust = -0.5, size = 3) +
  labs(
    title = "Distribución por categoría de edad del conductor",
    x = "Edad del conductor",
    y = "Número de pólizas"
  ) +
  expand_limits(y = max(porcentajes_edad$n) * 1.1) +  # aumenta el eje Y un 10%
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8)
  )

Tasa media de siniestros por tramo de edad:

importe_por_edad <- polizas %>%
  group_by(edad_conductor_cat) %>%
  summarise(importe_medio = mean(importe_medio_siniestro, na.rm = TRUE)) %>%
  ungroup()


ggplot(importe_por_edad,
       aes(x = edad_conductor_cat, y = importe_medio)) +
  geom_col(fill = "steelblue", color = "white") +
  geom_text(aes(label = sprintf("%.0f €", importe_medio)),
            vjust = -0.5, size = 3) +
  labs(
    title = "Importe medio del siniestro por grupo de edad del conductor",
    x = "Grupo de edad del conductor",
    y = "Importe medio (€)"
  ) +
  expand_limits(y = max(importe_por_edad$importe_medio) * 1.15) +  # ← aumenta eje Y un 15%
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.title = element_text(size = 10),
    axis.text  = element_text(size = 8),
    panel.grid.minor = element_blank()
  )

Edad del vehículo

Depuramos la base de datos:

# Modificar las observaciones con exposición > 1, igualándolas a 1
polizas <- polizas %>%
  mutate(edad_vehiculo = ifelse(edad_vehiculo > 30, 30, edad_vehiculo))

# Verificar el cambio
summary(polizas$edad_vehiculo)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   2.000   6.000   7.029  11.000  30.000

Distribución de la variable edad del vehículo depurada:

ggplot(polizas, aes(x = edad_vehiculo)) +
  geom_histogram(bins = 30, fill = "steelblue", color = "white") +
  labs(
    title = "Distribución de la edad del vehículo",
    x = "Años del vehículo",
    y = "Número de pólizas"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 12),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8)
  )

Categorizamos la variable

polizas <- polizas %>%
  mutate(
    edad_vehiculo_cat = case_when(
      edad_vehiculo <= 5  ~ "Nuevo",
      edad_vehiculo <= 12 ~ "Medio",
      TRUE                ~ "Viejo"
    )
  )

polizas$edad_vehiculo_cat <- factor(
  polizas$edad_vehiculo_cat,
  levels = c("Nuevo", "Medio", "Viejo")
)

Distribución variable edad del vehículo categorizada:

polizas %>%
  count(edad_vehiculo_cat) %>%
  mutate(porcentaje = 100 * n / sum(n)) %>%
  ggplot(aes(x = edad_vehiculo_cat, y = porcentaje)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = sprintf("%.1f%%", porcentaje)),
            vjust = -0.3, size = 3) +
  labs(
    title = "Distribución porcentual por categorías de edad del vehículo",
    x = "Categoría de edad del vehículo",
    y = "Porcentaje (%)"
  ) +
  theme_minimal(base_size = 10)

Tasa media de siniestros por categoría de edad del vehículo:

tasa_cat <- polizas %>%
  group_by(edad_vehiculo_cat) %>%
  summarise(
    tasa_media = mean(tasa_siniestros, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(tasa_cat, aes(x = edad_vehiculo_cat, y = tasa_media)) +
   geom_col(fill = "steelblue") +
  geom_text(aes(label = round(tasa_media, 3)),
            vjust = -0.3, size = 3) +
  labs(
    title = "Tasa media de siniestros según la edad del vehículo",
    x = "Categoría de edad del vehículo",
    y = "Tasa media de siniestros"
  )+
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 12),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8)
  )

Importe medio de siniestro por categoría del vehículo

importe_cat <- polizas %>%
  group_by(edad_vehiculo_cat) %>%
  summarise(
    importe_medio = mean(importe_medio_siniestro, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(importe_cat, aes(x = edad_vehiculo_cat, y = importe_medio)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = paste0(round(importe_medio, 0), " €")),
            vjust = -0.3, size = 3) +
  labs(
    title = "Importe medio de siniestro según la edad del vehículo",
    x = "Categoría de edad del vehículo",
    y = "Importe medio (€)"
  )+
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 12),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8))

Potencia

ggplot(polizas, aes(x = potencia)) +
  geom_bar(fill = "steelblue", color = "white") +
  labs(
    title = "Distribución de la potencia del vehículo",
    x = "Potencia",
    y = "Número de pólizas"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 12),
    axis.title = element_text(size = 10),
    axis.text  = element_text(size = 8)
  )

siniestros_medios_por_pot <- polizas %>%
  group_by(potencia) %>%
  summarise(
    siniestros_medios = mean(numero_siniestros, na.rm = TRUE)
  ) %>%
  ungroup()

ggplot(siniestros_medios_por_pot, aes(x = potencia, y = siniestros_medios)) +
  geom_col(fill = "steelblue", color = "white") +
  geom_text(aes(label = sprintf("%.3f", siniestros_medios)), 
            vjust = -0.5, size = 3) +
  labs(
    title = "Número medio de siniestros por potencia del vehículo",
    x = "Potencia del vehículo",
    y = "Número medio de siniestros por póliza"
  )+
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 12),
    axis.title = element_text(size = 10),
    axis.text  = element_text(size = 8)
  )

tasa_por_pot <- polizas %>%
  group_by(potencia) %>%
  summarise(
    tasa_media = sum(numero_siniestros, na.rm = TRUE) / sum(exposicion, na.rm = TRUE)
  ) %>%
  ungroup()

ggplot(tasa_por_pot, aes(x = potencia, y = tasa_media)) +
  geom_line(color = "steelblue", linewidth = 1.2) +
  geom_point(color = "steelblue", size = 2) +
  geom_text(aes(label = sprintf("%.3f", tasa_media)), 
            vjust = -0.8, size = 3) +
  labs(
    title = "Tasa media de siniestralidad por potencia del vehículo",
    x = "Potencia del vehículo",
    y = "Tasa media de siniestralidad"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 12),
    axis.title = element_text(size = 10),
    axis.text  = element_text(size = 8)
  )

Potencia en grupo

polizas <- polizas %>%
  mutate(
    potencia_grupo = case_when(
      potencia <= 5 ~ "Baja",
      potencia > 5 & potencia <= 9 ~ "Media",
      potencia > 9 & potencia <= 12 ~ "Alta",
      potencia > 12 ~ "Muy Alta",
      TRUE ~ NA_character_
    )
  )

ggplot(polizas, aes(x = potencia_grupo)) +
  geom_bar(fill = "steelblue", color = "white") +
  geom_text(stat = "count", aes(label = scales::comma(..count..)), 
            vjust = -0.4, size = 2.8) +   # añade etiquetas encima de las barras
  labs(
    title = "Distribución por grupo de potencia",
    x = "Grupo de potencia",
    y = "Número de pólizas"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 12),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8)
  )

## Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(count)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

### Modelización de la frecuencia

Preparamos y comprobamos la variable respuesta

summary(polizas$numero_siniestros)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.03895 0.00000 4.00000

table(polizas$numero_siniestros)

## 
##      0      1      2      3      4 
## 653047  23571   1298     62     13

mean(polizas$numero_siniestros)

## [1] 0.03894594

var(polizas$numero_siniestros)

## [1] 0.04203695

mean(polizas$numero_siniestros == 0)

## [1] 0.963209

Probamos la distribución Poisson

Función auxiliar para la sobredispersión:

dispersion <- function(model) {
  sum(residuals(model, type = "pearson")^2) / df.residual(model)
}

polizas$potencia_grupo <- factor(polizas$potencia_grupo, ordered = FALSE)

Ajustamos el modelo Poisson básico

modelo_posion_base <- glm(numero_siniestros ~ edad_conductor_cat + edad_vehiculo_cat + potencia_grupo + combustible + bonus_malus,
              offset = log(exposicion),
              family = poisson,
              data = polizas)

summary(modelo_posion_base)

## 
## Call:
## glm(formula = numero_siniestros ~ edad_conductor_cat + edad_vehiculo_cat + 
##     potencia_grupo + combustible + bonus_malus, family = poisson, 
##     data = polizas, offset = log(exposicion))
## 
## Coefficients:
##                                     Estimate Std. Error z value
## (Intercept)                       -4.1510498  0.0441119 -94.103
## edad_conductor_catAdultos jóvenes -0.0833998  0.0245830  -3.393
## edad_conductor_catAdultos          0.2175425  0.0256199   8.491
## edad_conductor_catMayores          0.1178927  0.0290862   4.053
## edad_vehiculo_catMedio             0.0521224  0.0136185   3.827
## edad_vehiculo_catViejo            -0.2121828  0.0177988 -11.921
## potencia_grupoBaja                -0.1932816  0.0238202  -8.114
## potencia_grupoMedia               -0.1261790  0.0227993  -5.534
## potencia_grupoMuy Alta             0.0859290  0.0586893   1.464
## combustibleRegular                -0.1313401  0.0126962 -10.345
## bonus_malus                        0.0274305  0.0003405  80.566
##                                               Pr(>|z|)    
## (Intercept)                       < 0.0000000000000002 ***
## edad_conductor_catAdultos jóvenes             0.000692 ***
## edad_conductor_catAdultos         < 0.0000000000000002 ***
## edad_conductor_catMayores         0.000050518233659135 ***
## edad_vehiculo_catMedio                        0.000130 ***
## edad_vehiculo_catViejo            < 0.0000000000000002 ***
## potencia_grupoBaja                0.000000000000000489 ***
## potencia_grupoMedia               0.000000031238680864 ***
## potencia_grupoMuy Alta                        0.143158    
## combustibleRegular                < 0.0000000000000002 ***
## bonus_malus                       < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 171311  on 677990  degrees of freedom
## Residual deviance: 164937  on 677980  degrees of freedom
## AIC: 215722
## 
## Number of Fisher Scoring iterations: 6

dispersion(modelo_posion_base)

## [1] 1.664483

AIC(modelo_posion_base)

## [1] 215721.7

par(mfrow = c(2, 2)); plot(modelo_posion_base); par(mfrow = c(1, 1))

Probamos la distribución Binomial negativa

modelo_binomial_negativo <- glm.nb(
  numero_siniestros ~ edad_conductor_cat +
                      edad_vehiculo_cat +
                      potencia_grupo +
                      combustible +
                      bonus_malus +
                      offset(log(exposicion)),
  data = polizas
)

summary(modelo_binomial_negativo)

## 
## Call:
## glm.nb(formula = numero_siniestros ~ edad_conductor_cat + edad_vehiculo_cat + 
##     potencia_grupo + combustible + bonus_malus + offset(log(exposicion)), 
##     data = polizas, init.theta = 1.043013931, link = log)
## 
## Coefficients:
##                                     Estimate Std. Error z value
## (Intercept)                       -4.1850273  0.0468518 -89.325
## edad_conductor_catAdultos jóvenes -0.0778844  0.0257736  -3.022
## edad_conductor_catAdultos          0.2290706  0.0270195   8.478
## edad_conductor_catMayores          0.1264174  0.0305558   4.137
## edad_vehiculo_catMedio             0.0539929  0.0140405   3.846
## edad_vehiculo_catViejo            -0.2128577  0.0182840 -11.642
## potencia_grupoBaja                -0.1907287  0.0245545  -7.768
## potencia_grupoMedia               -0.1248341  0.0235090  -5.310
## potencia_grupoMuy Alta             0.0892657  0.0604698   1.476
## combustibleRegular                -0.1320668  0.0130835 -10.094
## bonus_malus                        0.0278758  0.0003724  74.851
##                                               Pr(>|z|)    
## (Intercept)                       < 0.0000000000000002 ***
## edad_conductor_catAdultos jóvenes              0.00251 ** 
## edad_conductor_catAdultos         < 0.0000000000000002 ***
## edad_conductor_catMayores            0.000035146351850 ***
## edad_vehiculo_catMedio                         0.00012 ***
## edad_vehiculo_catViejo            < 0.0000000000000002 ***
## potencia_grupoBaja                   0.000000000000008 ***
## potencia_grupoMedia                  0.000000109589658 ***
## potencia_grupoMuy Alta                         0.13989    
## combustibleRegular                < 0.0000000000000002 ***
## bonus_malus                       < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.043) family taken to be 1)
## 
##     Null deviance: 151138  on 677990  degrees of freedom
## Residual deviance: 145155  on 677980  degrees of freedom
## AIC: 215185
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.0430 
##           Std. Err.:  0.0576 
## 
##  2 x log-likelihood:  -215161.2840

AIC(modelo_binomial_negativo)

## [1] 215185.3

dispersion(modelo_binomial_negativo)

## [1] 1.622827

par(mfrow = c(2, 2)); plot(modelo_binomial_negativo); par(mfrow = c(1, 1))

Estructura del predictor usando Poisson

Edad del conductor

1. Modelo poisson con edad Continua lineal

m_pois_age_lin <- glm(
  numero_siniestros ~ edad_conductor +
                      edad_vehiculo_cat +
                      potencia_grupo +
                      combustible +
                      bonus_malus +
                      offset(log(exposicion)),
  family = poisson,
  data   = polizas
)

AIC(modelo_posion_base, m_pois_age_lin)

##                    df      AIC
## modelo_posion_base 11 215721.7
## m_pois_age_lin      9 215985.6

dispersion(m_pois_age_lin)

## [1] 1.67128

2. Edad cuadrática

m_pois_age_quad <- glm(
  numero_siniestros ~ edad_conductor + I(edad_conductor^2) +
                      edad_vehiculo_cat +
                      potencia_grupo +
                      combustible +
                      bonus_malus +
                      offset(log(exposicion)),
  family = poisson,
  data   = polizas
)

AIC(modelo_posion_base, m_pois_age_quad)

##                    df      AIC
## modelo_posion_base 11 215721.7
## m_pois_age_quad    10 215890.0

dispersion(m_pois_age_quad)

## [1] 1.665208

3. Edad con spline

library(splines)

m_pois_age_spline <- glm(
  numero_siniestros ~ ns(edad_conductor, df = 4) +
                      edad_vehiculo_cat +
                      potencia_grupo +
                      combustible +
                      bonus_malus +
                      offset(log(exposicion)),
  family = poisson,
  data   = polizas
)

AIC(modelo_posion_base, m_pois_age_spline)

##                    df      AIC
## modelo_posion_base 11 215721.7
## m_pois_age_spline  12 215544.5

dispersion(m_pois_age_spline)

## [1] 1.661709

Edad vehículo

1. Edad vehículo continua

m_pois_veh_lin <- glm(
  numero_siniestros ~ edad_conductor_cat +   # o tu forma elegida arriba
                      edad_vehiculo +
                      potencia_grupo +
                      combustible +
                      bonus_malus +
                      offset(log(exposicion)),
  family = poisson,
  data   = polizas
)

AIC(m_pois_veh_lin)

## [1] 215840.3

dispersion(m_pois_veh_lin)

## [1] 1.659176

2. Edad vehículo cuadrática

m_pois_veh_quad <- glm(
  numero_siniestros ~ edad_conductor_cat +
                      edad_vehiculo + I(edad_vehiculo^2) +
                      potencia_grupo +
                      combustible +
                      bonus_malus +
                      offset(log(exposicion)),
  family = poisson,
  data = polizas
)

AIC(m_pois_veh_quad)

## [1] 215676.7

dispersion(m_pois_veh_quad)

## [1] 1.662875

m_pois_veh_spline <- glm(
  numero_siniestros ~ edad_conductor_cat +
                      ns(edad_vehiculo, df = 4) +
                      potencia_grupo +
                      combustible +
                      bonus_malus +
                      offset(log(exposicion)),
  family = poisson,
  data = polizas
)

AIC(m_pois_veh_spline)

## [1] 215661

dispersion(m_pois_veh_spline)

## [1] 1.66449

Potencia

1. Potencia lineal

m_pois_pot_lin <- glm(
  numero_siniestros ~ edad_conductor_cat +
                      edad_vehiculo_cat +
                      potencia +
                      combustible +
                      bonus_malus +
                      offset(log(exposicion)),
  family = poisson,
  data = polizas
)

AIC(m_pois_pot_lin)

## [1] 215696.6

dispersion(m_pois_pot_lin)

## [1] 1.665391

m_pois_veh_quad_lin <- glm(
  numero_siniestros ~ edad_conductor_cat +
                      edad_vehiculo + I(edad_vehiculo^2) +
                      potencia +
                      combustible +
                      bonus_malus +
                      offset(log(exposicion)),
  family = poisson,
  data = polizas
)

AIC(m_pois_veh_quad_lin)

## [1] 215648.9

dispersion(m_pois_veh_quad_lin)

## [1] 1.664081

2. Potencia cuadrática

m_pois_pot_quad <- glm(
  numero_siniestros ~ edad_conductor_cat +
                      edad_vehiculo_cat +
                      potencia + I(potencia^2) +
                      combustible +
                      bonus_malus +
                      offset(log(exposicion)),
  family = poisson,
  data = polizas
)

AIC(m_pois_pot_quad)

## [1] 215698.1

dispersion(m_pois_pot_quad)

## [1] 1.665872

3. Potencia Log

m_pois_pot_log <- glm(
  numero_siniestros ~ edad_conductor_cat +
                      edad_vehiculo_cat +
                      log(potencia) +
                      combustible +
                      bonus_malus +
                      offset(log(exposicion)),
  family = poisson,
  data = polizas
)

AIC(m_pois_pot_log)

## [1] 215696.9

dispersion(m_pois_pot_log)

## [1] 1.666867

Actualizo modelo base

modelo_poison_base_1 <- glm(
  numero_siniestros ~ edad_conductor_cat +
                      edad_vehiculo + I(edad_vehiculo^2) +
                      potencia +
                      combustible +
                      bonus_malus +
                      offset(log(exposicion)),
  family = poisson,
  data = polizas
)

summary(modelo_poison_base_1)

## 
## Call:
## glm(formula = numero_siniestros ~ edad_conductor_cat + edad_vehiculo + 
##     I(edad_vehiculo^2) + potencia + combustible + bonus_malus + 
##     offset(log(exposicion)), family = poisson, data = polizas)
## 
## Coefficients:
##                                     Estimate Std. Error  z value
## (Intercept)                       -4.5529118  0.0454218 -100.236
## edad_conductor_catAdultos jóvenes -0.0852246  0.0245875   -3.466
## edad_conductor_catAdultos          0.2150798  0.0256282    8.392
## edad_conductor_catMayores          0.1201094  0.0290828    4.130
## edad_vehiculo                      0.0323555  0.0036902    8.768
## I(edad_vehiculo^2)                -0.0025241  0.0002015  -12.529
## potencia                           0.0338055  0.0030717   11.005
## combustibleRegular                -0.1336424  0.0124918  -10.698
## bonus_malus                        0.0274244  0.0003404   80.567
##                                               Pr(>|z|)    
## (Intercept)                       < 0.0000000000000002 ***
## edad_conductor_catAdultos jóvenes             0.000528 ***
## edad_conductor_catAdultos         < 0.0000000000000002 ***
## edad_conductor_catMayores                    0.0000363 ***
## edad_vehiculo                     < 0.0000000000000002 ***
## I(edad_vehiculo^2)                < 0.0000000000000002 ***
## potencia                          < 0.0000000000000002 ***
## combustibleRegular                < 0.0000000000000002 ***
## bonus_malus                       < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 171311  on 677990  degrees of freedom
## Residual deviance: 164868  on 677982  degrees of freedom
## AIC: 215649
## 
## Number of Fisher Scoring iterations: 6

AIC(modelo_poison_base_1)

## [1] 215648.9

dispersion(modelo_poison_base_1)

## [1] 1.664081

Combustible

1. Estudio si se queda o no la variable

summary(modelo_poison_base_1)$coefficients["combustibleRegular",]

##                             Estimate                           Std. Error 
##  -0.13364239045306164355153555334255   0.01249178909848610782851174860753 
##                              z value                             Pr(>|z|) 
## -10.69841872924815007195320504251868   0.00000000000000000000000001035298

m_pois_no_comb <- glm(
  numero_siniestros ~ edad_conductor_cat +
                      edad_vehiculo_cat +
                      potencia +
                      bonus_malus +
                      offset(log(exposicion)),
  family = poisson,
  data = polizas
)

AIC(modelo_poison_base_1, m_pois_no_comb)

##                      df      AIC
## modelo_poison_base_1  9 215648.9
## m_pois_no_comb        8 215808.6

dispersion(m_pois_no_comb)

## [1] 1.65824

anova(modelo_poison_base_1, m_pois_no_comb, test="Chisq")

## Analysis of Deviance Table
## 
## Model 1: numero_siniestros ~ edad_conductor_cat + edad_vehiculo + I(edad_vehiculo^2) + 
##     potencia + combustible + bonus_malus + offset(log(exposicion))
## Model 2: numero_siniestros ~ edad_conductor_cat + edad_vehiculo_cat + 
##     potencia + bonus_malus + offset(log(exposicion))
##   Resid. Df Resid. Dev Df Deviance              Pr(>Chi)    
## 1    677982     164868                                      
## 2    677983     165030 -1  -161.71 < 0.00000000000000022 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2. Pruebo interacción

m_pois_interact <- glm(
  numero_siniestros ~ edad_conductor_cat +
                       edad_vehiculo + I(edad_vehiculo^2)+
                      potencia*combustible +
                      bonus_malus +
                      offset(log(exposicion)),
  family = poisson,
  data   = polizas
)

AIC(modelo_poison_base_1, m_pois_interact)

##                      df      AIC
## modelo_poison_base_1  9 215648.9
## m_pois_interact      10 215635.3

dispersion(m_pois_interact)

## [1] 1.66635

anova(modelo_poison_base_1, m_pois_interact, test="Chisq")

## Analysis of Deviance Table
## 
## Model 1: numero_siniestros ~ edad_conductor_cat + edad_vehiculo + I(edad_vehiculo^2) + 
##     potencia + combustible + bonus_malus + offset(log(exposicion))
## Model 2: numero_siniestros ~ edad_conductor_cat + edad_vehiculo + I(edad_vehiculo^2) + 
##     potencia * combustible + bonus_malus + offset(log(exposicion))
##   Resid. Df Resid. Dev Df Deviance   Pr(>Chi)    
## 1    677982     164868                           
## 2    677981     164853  1   15.562 0.00007983 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Actualizo modelo base

m_pois_final <-  glm(
  numero_siniestros ~ edad_conductor_cat +
                       edad_vehiculo + I(edad_vehiculo^2)+
                      potencia*combustible +
                      bonus_malus +
                      offset(log(exposicion)),
  family = poisson,
  data   = polizas
)

summary(m_pois_final)

## 
## Call:
## glm(formula = numero_siniestros ~ edad_conductor_cat + edad_vehiculo + 
##     I(edad_vehiculo^2) + potencia * combustible + bonus_malus + 
##     offset(log(exposicion)), family = poisson, data = polizas)
## 
## Coefficients:
##                                     Estimate Std. Error z value
## (Intercept)                       -4.4547759  0.0518204 -85.966
## edad_conductor_catAdultos jóvenes -0.0839910  0.0245884  -3.416
## edad_conductor_catAdultos          0.2192898  0.0256506   8.549
## edad_conductor_catMayores          0.1237228  0.0290992   4.252
## edad_vehiculo                      0.0313359  0.0036940   8.483
## I(edad_vehiculo^2)                -0.0024707  0.0002016 -12.258
## potencia                           0.0189153  0.0049032   3.858
## combustibleRegular                -0.2935956  0.0424651  -6.914
## bonus_malus                        0.0274189  0.0003404  80.543
## potencia:combustibleRegular        0.0246730  0.0062642   3.939
##                                               Pr(>|z|)    
## (Intercept)                       < 0.0000000000000002 ***
## edad_conductor_catAdultos jóvenes             0.000636 ***
## edad_conductor_catAdultos         < 0.0000000000000002 ***
## edad_conductor_catMayores             0.00002120980838 ***
## edad_vehiculo                     < 0.0000000000000002 ***
## I(edad_vehiculo^2)                < 0.0000000000000002 ***
## potencia                                      0.000114 ***
## combustibleRegular                    0.00000000000472 ***
## bonus_malus                       < 0.0000000000000002 ***
## potencia:combustibleRegular           0.00008190507532 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 171311  on 677990  degrees of freedom
## Residual deviance: 164853  on 677981  degrees of freedom
## AIC: 215635
## 
## Number of Fisher Scoring iterations: 6

AIC(m_pois_final)

## [1] 215635.3

dispersion(m_pois_final)

## [1] 1.66635

par(mfrow=c(2,2)); plot(m_pois_final); par(mfrow=c(1,1))

Ajustamos Binomial negativa

m_nb_final <- glm.nb(
  numero_siniestros ~ edad_conductor_cat +
                       edad_vehiculo + I(edad_vehiculo^2)+
                      potencia*combustible +
                      bonus_malus +
                      offset(log(exposicion)),
  data = polizas
)

summary(m_nb_final)

## 
## Call:
## glm.nb(formula = numero_siniestros ~ edad_conductor_cat + edad_vehiculo + 
##     I(edad_vehiculo^2) + potencia * combustible + bonus_malus + 
##     offset(log(exposicion)), data = polizas, init.theta = 1.048476778, 
##     link = log)
## 
## Coefficients:
##                                     Estimate Std. Error z value
## (Intercept)                       -4.4838603  0.0546565 -82.037
## edad_conductor_catAdultos jóvenes -0.0783767  0.0257741  -3.041
## edad_conductor_catAdultos          0.2307577  0.0270422   8.533
## edad_conductor_catMayores          0.1325212  0.0305647   4.336
## edad_vehiculo                      0.0315287  0.0037834   8.334
## I(edad_vehiculo^2)                -0.0024803  0.0002061 -12.036
## potencia                           0.0184044  0.0050530   3.642
## combustibleRegular                -0.2969096  0.0437502  -6.786
## bonus_malus                        0.0278638  0.0003723  74.845
## potencia:combustibleRegular        0.0250866  0.0064541   3.887
##                                               Pr(>|z|)    
## (Intercept)                       < 0.0000000000000002 ***
## edad_conductor_catAdultos jóvenes             0.002359 ** 
## edad_conductor_catAdultos         < 0.0000000000000002 ***
## edad_conductor_catMayores              0.0000145261777 ***
## edad_vehiculo                     < 0.0000000000000002 ***
## I(edad_vehiculo^2)                < 0.0000000000000002 ***
## potencia                                      0.000270 ***
## combustibleRegular                     0.0000000000115 ***
## bonus_malus                       < 0.0000000000000002 ***
## potencia:combustibleRegular                   0.000102 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.0485) family taken to be 1)
## 
##     Null deviance: 151219  on 677990  degrees of freedom
## Residual deviance: 145153  on 677981  degrees of freedom
## AIC: 215102
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.0485 
##           Std. Err.:  0.0581 
## 
##  2 x log-likelihood:  -215080.3840

comparamos con Poisson

tabla1 <- AIC(m_pois_final, m_nb_final)
dispersion(m_nb_final)

## [1] 1.624849

Revisamos los residuos

par(mfrow=c(2,2)); plot(m_nb_final); par(mfrow=c(1,1))

Diagnóstico del modelo NB

Gráficos básicos de diagnóstico

par(mfrow = c(2, 2))
plot(m_nb_final)    

par(mfrow = c(1, 1))

res_nb  <- residuals(m_nb_final, type = "pearson")
lev_nb  <- hatvalues(m_nb_final)
cook_nb <- cooks.distance(m_nb_final)

sort(cook_nb, decreasing=TRUE)[1:10]

##       259414        39811       141089       558610       574655       152687 
## 0.0032881621 0.0031219588 0.0025989982 0.0014153098 0.0009295472 0.0009032286 
##       469045       541664       589135       259341 
## 0.0009012370 0.0008936455 0.0008918052 0.0008802308

eff_df <- as.data.frame(effect("edad_vehiculo", m_nb_final))

ggplot(eff_df, aes(x = edad_vehiculo, y = fit)) +
  geom_line(fill = "steelblue",size = 1) +
  geom_ribbon(fill = "steelblue",aes(ymin = lower, ymax = upper), alpha = 0.2) +
  labs(
    title = "Efecto no lineal de la edad del vehículo sobre la frecuencia siniestral",
    x = "Edad del vehículo",
    y = "Frecuencia esperada de siniestros"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 12),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8)
  )

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning in geom_line(fill = "steelblue", size = 1): Ignoring unknown
## parameters: `fill`

Modelización de la severidad

Definir la variable severidad y conjunto de datos

datos_sev <- polizas %>%
  filter(numero_siniestros > 0) %>%
  mutate(
    sev_media = total_importe_siniestros / numero_siniestros
  )

summary(datos_sev$sev_media)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##       1.0     710.6    1172.0    2222.3    1228.1 4075400.6

quantile(datos_sev$sev_media, probs = c(0.90, 0.95, 0.99, 0.995, 0.999), na.rm = TRUE)

##        90%        95%        99%      99.5%      99.9% 
##   2742.509   4659.463  16327.001  34564.646 133304.840

corte <- quantile(datos_sev$sev_media, probs = 0.995, na.rm = TRUE)

datos_sev$sev_trunc <- pmin(datos_sev$sev_media, corte)

summary(datos_sev$sev_trunc)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     1.0   710.6  1172.0  1683.8  1228.1 34564.6

prop_trunc <- mean(datos_sev$sev_media > corte)
prop_trunc

## [1] 0.005011225

Exploratorio rápido de la severidad (para justificar Gamma/log)

ggplot(datos_sev, aes(x = sev_trunc)) +
  geom_histogram(bins = 50, fill = "#4C72B0", alpha = 0.9) +
  labs(
    title = "Severidad truncada",
    x = "Severidad",
    y = "Número de pólizas"
  ) +  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 12),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8)
  )

ggplot(datos_sev, aes(x = log(sev_trunc))) +
  geom_histogram(bins = 50, fill = "#4C72B0", alpha = 0.9) +
  labs(
    title = "Log-severidad",
    x = "log(Severidad)",
    y = "Número de pólizas"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 12),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8)
  )

Definimos los factores explicativos

datos_sev <- datos_sev %>%
  mutate(
    edad_conductor_cat = factor(edad_conductor_cat),
    edad_vehiculo_cat  = factor(edad_vehiculo_cat),
    potencia_grupo     = factor(potencia_grupo),
    combustible        = factor(combustible),
    bonus_malus_num = as.numeric(as.character(bonus_malus))
  )

Selección del predictor lineal

mod_sev_0 <- glm(
  sev_trunc ~ edad_conductor_cat +
              edad_vehiculo_cat +
              potencia_grupo +
              bonus_malus,
  family  = Gamma(link = "log"),
  data    = datos_sev,
  weights = numero_siniestros
)
summary(mod_sev_0)

## 
## Call:
## glm(formula = sev_trunc ~ edad_conductor_cat + edad_vehiculo_cat + 
##     potencia_grupo + bonus_malus, family = Gamma(link = "log"), 
##     data = datos_sev, weights = numero_siniestros)
## 
## Coefficients:
##                                     Estimate Std. Error t value
## (Intercept)                        7.4775042  0.0900126  83.072
## edad_conductor_catAdultos jóvenes -0.0846146  0.0490522  -1.725
## edad_conductor_catAdultos         -0.0732701  0.0516660  -1.418
## edad_conductor_catMayores         -0.0292565  0.0578605  -0.506
## edad_vehiculo_catMedio            -0.0901533  0.0266370  -3.385
## edad_vehiculo_catViejo            -0.1126957  0.0346615  -3.251
## potencia_grupoBaja                -0.1150261  0.0462989  -2.484
## potencia_grupoMedia               -0.0502134  0.0445056  -1.128
## potencia_grupoMuy Alta             0.0677571  0.1141975   0.593
## bonus_malus                        0.0023140  0.0007252   3.191
##                                               Pr(>|t|)    
## (Intercept)                       < 0.0000000000000002 ***
## edad_conductor_catAdultos jóvenes             0.084541 .  
## edad_conductor_catAdultos                     0.156159    
## edad_conductor_catMayores                     0.613115    
## edad_vehiculo_catMedio                        0.000714 ***
## edad_vehiculo_catViejo                        0.001150 ** 
## potencia_grupoBaja                            0.012983 *  
## potencia_grupoMedia                           0.259225    
## potencia_grupoMuy Alta                        0.552964    
## bonus_malus                                   0.001420 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 3.806349)
## 
##     Null deviance: 30402  on 24943  degrees of freedom
## Residual deviance: 30235  on 24934  degrees of freedom
## AIC: 446180
## 
## Number of Fisher Scoring iterations: 6

Añadimos combustible

mod_sev_1 <- glm(
  sev_trunc ~ edad_conductor_cat +
              edad_vehiculo_cat +
              potencia_grupo +
              bonus_malus +
              combustible,
  family  = Gamma(link = "log"),
  data    = datos_sev,
  weights = numero_siniestros
)
summary(mod_sev_1)

## 
## Call:
## glm(formula = sev_trunc ~ edad_conductor_cat + edad_vehiculo_cat + 
##     potencia_grupo + bonus_malus + combustible, family = Gamma(link = "log"), 
##     data = datos_sev, weights = numero_siniestros)
## 
## Coefficients:
##                                     Estimate Std. Error t value
## (Intercept)                        7.4696116  0.0904991  82.538
## edad_conductor_catAdultos jóvenes -0.0833713  0.0490685  -1.699
## edad_conductor_catAdultos         -0.0727989  0.0516524  -1.409
## edad_conductor_catMayores         -0.0309608  0.0578777  -0.535
## edad_vehiculo_catMedio            -0.0910682  0.0266443  -3.418
## edad_vehiculo_catViejo            -0.1164238  0.0349163  -3.334
## potencia_grupoBaja                -0.1166356  0.0463453  -2.517
## potencia_grupoMedia               -0.0489151  0.0445077  -1.099
## potencia_grupoMuy Alta             0.0599634  0.1144862   0.524
## bonus_malus                        0.0022927  0.0007253   3.161
## combustibleRegular                 0.0200524  0.0246449   0.814
##                                               Pr(>|t|)    
## (Intercept)                       < 0.0000000000000002 ***
## edad_conductor_catAdultos jóvenes             0.089316 .  
## edad_conductor_catAdultos                     0.158729    
## edad_conductor_catMayores                     0.592700    
## edad_vehiculo_catMedio                        0.000632 ***
## edad_vehiculo_catViejo                        0.000856 ***
## potencia_grupoBaja                            0.011853 *  
## potencia_grupoMedia                           0.271768    
## potencia_grupoMuy Alta                        0.600450    
## bonus_malus                                   0.001575 ** 
## combustibleRegular                            0.415853    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 3.803772)
## 
##     Null deviance: 30402  on 24943  degrees of freedom
## Residual deviance: 30232  on 24933  degrees of freedom
## AIC: 446179
## 
## Number of Fisher Scoring iterations: 6

AIC(mod_sev_0, mod_sev_1)

##           df      AIC
## mod_sev_0 11 446179.9
## mod_sev_1 12 446179.3

Modelo final

mod_sev_final <- mod_sev_0  

Diagnóstico del modelo final Revisamos la dispersión y calidad global del modelo

dispersion <- mod_sev_final$deviance / mod_sev_final$df.residual
dispersion

## [1] 1.212588

Interpretar los coeficientes

coef_sev <- coef(mod_sev_final)
relatividades <- exp(coef_sev)
round(relatividades, 3)

##                       (Intercept) edad_conductor_catAdultos jóvenes 
##                          1767.823                             0.919 
##         edad_conductor_catAdultos         edad_conductor_catMayores 
##                             0.929                             0.971 
##            edad_vehiculo_catMedio            edad_vehiculo_catViejo 
##                             0.914                             0.893 
##                potencia_grupoBaja               potencia_grupoMedia 
##                             0.891                             0.951 
##            potencia_grupoMuy Alta                       bonus_malus 
##                             1.070                             1.002

Diagnóstico del modelo de severidad

# Preparar datos de diagnóstico
datos_diag <- datos_sev %>%
  mutate(
    fitted = fitted(mod_sev_0),
    resid_dev = residuals(mod_sev_0, type = "deviance"),
    resid_pearson = residuals(mod_sev_0, type = "pearson"),
    cook = cooks.distance(mod_sev_0)
  )

# Residuos vs Ajustados
p1 <- ggplot(datos_diag, aes(x = fitted, y = resid_dev)) +
  geom_point(alpha = 0.25) +
  geom_hline(yintercept = 0, linetype = 2) +
  geom_smooth(se = FALSE, method = "loess") +
  labs(title = "Residuos devianza vs ajustados", x = "μ̂", y = "Residuos (dev)") +
  theme_minimal()

# Q-Q plot residuos Pearson
p2 <- ggplot(datos_diag %>% filter(is.finite(resid_pearson)),
             aes(sample = resid_pearson)) +
  stat_qq(alpha = 0.25) +
  stat_qq_line() +
  labs(title = "Q-Q residuos Pearson", x = "Cuantiles teóricos", y = "Residuos") +
  theme_minimal()

# Cook's distance
n <- nobs(mod_sev_0)
umbral <- 4 / n
p3 <- ggplot(datos_diag, aes(x = seq_along(cook), y = cook)) +
  geom_col() +
  geom_hline(yintercept = umbral, linetype = 2) +
  labs(title = "Distancia de Cook", x = "Observación", y = "Cook") +
  theme_minimal()

# Calibración por deciles (O/E)
datos_cal <- datos_sev %>%
  mutate(
    sev_esp = fitted(mod_sev_0),
    sev_obs = sev_trunc,
    w = numero_siniestros
  ) %>%
  filter(is.finite(sev_esp), is.finite(sev_obs), w > 0) %>%
  mutate(decil = ntile(sev_esp, 10)) %>%
  group_by(decil) %>%
  summarise(
    sev_obs = sum(sev_obs * w) / sum(w),
    sev_esp = sum(sev_esp * w) / sum(w),
    OE = sev_obs / sev_esp,
    .groups = "drop"
  )

p4 <- ggplot(datos_cal, aes(x = decil, y = OE)) +
  geom_line() +
  geom_point(size = 2) +
  geom_hline(yintercept = 1, linetype = 2) +
  scale_x_continuous(breaks = 1:10) +
  labs(title = "Calibración por deciles (O/E)", x = "Decil μ̂", y = "O/E") +
  theme_minimal()

Modelo de Tweedie

Asiganamos la variable respuesta

polizas <- polizas %>%
  mutate(
    S_bruto = ifelse(numero_siniestros == 0, 0,
                     total_importe_siniestros)  
  )

summary(polizas$S_bruto)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##       0.0       0.0       0.0      88.4       0.0 4075400.6

mean(polizas$S_bruto == 0)

## [1] 0.963209

p_trunc <- 0.995

c_trunc <- quantile(
  polizas$S_bruto[polizas$S_bruto > 0],
  probs = p_trunc,
  na.rm = TRUE
)

c_trunc

##    99.5% 
## 35630.33

polizas <- polizas %>%
  mutate(
    S_trunc = pmin(S_bruto, as.numeric(c_trunc))
  )

summary(polizas$S_trunc)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##     0.00     0.00     0.00    66.13     0.00 35630.33

max(polizas$S_trunc)

## [1] 35630.33

c(
  umbral_truncado = as.numeric(c_trunc),
  proporcion_truncadas = mean(polizas$S_bruto > c_trunc),
  impacto_medio = mean(polizas$S_bruto - polizas$S_trunc)
)

##      umbral_truncado proporcion_truncadas        impacto_medio 
##     35630.3308000000         0.0001843682        22.2357667581

** Especificación y estimación del modelo Tweedie**

m_tw <- gam(
  S_trunc ~ edad_conductor_cat +
            edad_vehiculo_cat +
            potencia_grupo +
            bonus_malus +
            combustible +
            offset(log(exposicion)),
  family = tw(link = "log"),
  data   = polizas,
  method = "REML"
)

summary(m_tw)

## 
## Family: Tweedie(p=1.496) 
## Link function: log 
## 
## Formula:
## S_trunc ~ edad_conductor_cat + edad_vehiculo_cat + potencia_grupo + 
##     bonus_malus + combustible + offset(log(exposicion))
## 
## Parametric coefficients:
##                                     Estimate Std. Error t value
## (Intercept)                        3.5630992  0.0696660  51.145
## edad_conductor_catAdultos jóvenes -0.2791880  0.0376075  -7.424
## edad_conductor_catAdultos         -0.0021799  0.0398750  -0.055
## edad_conductor_catMayores         -0.0786176  0.0444586  -1.768
## edad_vehiculo_catMedio             0.0345868  0.0197314   1.753
## edad_vehiculo_catViejo            -0.2490611  0.0249864  -9.968
## potencia_grupoBaja                -0.2110242  0.0343368  -6.146
## potencia_grupoMedia               -0.0920358  0.0329181  -2.796
## potencia_grupoMuy Alta             0.0818471  0.0851122   0.962
## bonus_malus                        0.0290061  0.0005877  49.356
## combustibleRegular                -0.0847245  0.0182546  -4.641
##                                               Pr(>|t|)    
## (Intercept)                       < 0.0000000000000002 ***
## edad_conductor_catAdultos jóvenes    0.000000000000114 ***
## edad_conductor_catAdultos                      0.95640    
## edad_conductor_catMayores                      0.07701 .  
## edad_vehiculo_catMedio                         0.07962 .  
## edad_vehiculo_catViejo            < 0.0000000000000002 ***
## potencia_grupoBaja                   0.000000000796443 ***
## potencia_grupoMedia                            0.00518 ** 
## potencia_grupoMuy Alta                         0.33623    
## bonus_malus                       < 0.0000000000000002 ***
## combustibleRegular                   0.000003463380605 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## R-sq.(adj) =  0.000806   Deviance explained = 3.88%
## -REML = 3.1899e+05  Scale est. = 439.36    n = 677991

Diagnóstico global del ajuste

par(mfrow = c(2,2))
gam.check(m_tw)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 8 iterations.
## Gradient range [-0.08276832,-0.07467906]
## (score 318991.9 & scale 439.3561).
## Hessian positive definite, eigenvalue range [11029.18,66471.63].
## Model rank =  11 / 11

# Extra: residuos vs ajustados (deviance residuals)
fitted_mu <- predict(m_tw, type = "response")
res_dev   <- residuals(m_tw, type = "deviance")

plot(fitted_mu, res_dev,
     xlab = "Valores ajustados (μ̂)",
     ylab = "Residuos de devianza",
     main = "Residuos vs Ajustados")
abline(h = 0, lty = 2)

Estimación de la prima pura

frecuencia x severidad

polizas$mu_freq <- predict(m_nb_final, newdata = polizas, type = "response")

polizas$mu_sev <- predict(mod_sev_final, newdata = polizas, type = "response")

polizas$pp_fs <- polizas$mu_freq * polizas$mu_sev

summary(polizas$pp_fs)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##    0.0562   20.8027   57.3075   67.0712   93.4853 6923.8416

quantile(polizas$pp_fs, probs = c(0.5, 0.75, 0.9, 0.95, 0.99))

##       50%       75%       90%       95%       99% 
##  57.30750  93.48527 123.45410 165.34326 300.96726

library(ggplot2)

ggplot(polizas, aes(x = log(pp_fs))) +
  geom_histogram(bins = 50, fill = "steelblue", color = "white") +
  labs(
    title = "Distribución de la prima pura estimada (Frecuencia–Severidad)",
    x = "log(Prima pura estimada)",
    y = "Frecuencia"
  )  +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 12),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8)
  )

df_conc <- polizas %>%
  filter(is.finite(pp_fs)) %>%
  arrange(desc(pp_fs)) %>%
  mutate(
    prop_polizas = row_number() / n(),
    prop_coste = cumsum(pp_fs) / sum(pp_fs)
  )

ggplot(df_conc, aes(x = prop_polizas, y = prop_coste)) +
  geom_line(linewidth = 1, color = "black") +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "steelblue") +
  labs(
    title = "Concentración del coste esperado",
    x = "Proporción acumulada de pólizas (ordenadas de mayor a menor)",
    y = "Proporción acumulada del coste esperado"
  )  +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 12),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8)
  )

modelo compuesto

polizas <- droplevels(polizas)

mod_tw_final <- gam(
  S_trunc ~ edad_conductor_cat +
            edad_vehiculo_cat +
            potencia_grupo +
            bonus_malus +
            combustible +
            offset(log(exposicion)),
  family = tw(link = "log"),
  data   = polizas,
  method = "REML"
)

polizas$pp_tw <- predict(mod_tw_final, newdata = polizas, type = "response")

summary(polizas$pp_tw)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##    0.1803   25.6107   70.8986   82.1729  114.9888 6148.0818

ggplot(polizas, aes(x = log(pp_tw))) +
  geom_histogram(bins = 50, fill = "steelblue", color = "white") +
  labs(
    title = "Distribución de la prima pura estimada (Tweedie)",
    x = "log(Prima pura estimada)",
    y = "Frecuencia"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 12),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 8)
  )

Comparaciones

mean(polizas$pp_fs)

## [1] 67.07117

mean(polizas$pp_tw)

## [1] 82.1729

sum(polizas$pp_fs)

## [1] 45473647

sum(polizas$pp_tw)

## [1] 55712485

tabla_cuantiles <- rbind(
  FS = quantile(polizas$pp_fs, probs = c(0.5, 0.75, 0.9, 0.95, 0.99)),
  TW = quantile(polizas$pp_tw, probs = c(0.5, 0.75, 0.9, 0.95, 0.99))
)

tabla_cuantiles

##         50%       75%      90%      95%      99%
## FS 57.30750  93.48527 123.4541 165.3433 300.9673
## TW 70.89863 114.98877 146.2145 203.3443 373.2328

cor(polizas$pp_fs, polizas$pp_tw)

## [1] 0.9874853

polizas$pp_emp <- polizas$importe_medio_siniestro / polizas$exposicion


comp_agregada <- polizas %>%
  summarise(
    n = n(),
    O_total = sum(pp_emp, na.rm = TRUE),
    E_FS_total = sum(pp_fs, na.rm = TRUE),
    E_TW_total = sum(pp_tw, na.rm = TRUE),
    O_media = mean(pp_emp, na.rm = TRUE),
    E_FS_media = mean(pp_fs, na.rm = TRUE),
    E_TW_media = mean(pp_tw, na.rm = TRUE),
    O_mediana = median(pp_emp, na.rm = TRUE),
    E_FS_mediana = median(pp_fs, na.rm = TRUE),
    E_TW_mediana = median(pp_tw, na.rm = TRUE)
  )

comp_agregada

##        n   O_total E_FS_total E_TW_total O_media E_FS_media E_TW_media
## 1 677991 243491667   45473647   55712485 359.137   67.07117    82.1729
##   O_mediana E_FS_mediana E_TW_mediana
## 1         0      57.3075     70.89863

cal_fs <- polizas %>%
  filter(is.finite(pp_fs), is.finite(pp_emp)) %>%
  mutate(decil = ntile(pp_fs, 10)) %>%
  group_by(decil) %>%
  summarise(
    n = n(),
    O = mean(pp_emp, na.rm = TRUE),
    E = mean(pp_fs, na.rm = TRUE),
    OE = O / E
  )

cal_tw <- polizas %>%
  filter(is.finite(pp_tw), is.finite(pp_emp)) %>%
  mutate(decil = ntile(pp_tw, 10)) %>%
  group_by(decil) %>%
  summarise(
    n = n(),
    O = mean(pp_emp, na.rm = TRUE),
    E = mean(pp_tw, na.rm = TRUE),
    OE = O / E
  )

ggplot(cal_fs, aes(x = decil, y = OE)) +
  geom_line() + geom_point() +
  geom_hline(yintercept = 1, linetype = "dashed") +
  labs(title = "Calibración por deciles (FS)", x = "Decil (PP_FS)", y = "O/E") +
  theme_minimal()

ggplot(cal_tw, aes(x = decil, y = OE)) +
  geom_line() + geom_point() +
  geom_hline(yintercept = 1, linetype = "dashed") +
  labs(title = "Calibración por deciles (Tweedie)", x = "Decil (PP_TW)", y = "O/E") +
  theme_minimal()

Exportación datos dashboard

write.xlsx(polizas, "base_datos.xlsx")